###################################################
# Replication script for:                         #  
# "Policy Monitoring and Ministerial Survival:    #
# Evidence from a Multiparty Presidentialism"     #
# Authors: Thiago N. Silva and Alejandro Medina   #
###################################################

### Installing packages 
# install.packages(c("survival", "survminer", "dplyr","ggplot2", "ggthemes", "ggfortify", "gridExtra", "stargazer", "glue","frailtypack"))

### Loading packages
library(survival)
library(survminer)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(ggfortify)
library(gridExtra)
library(stargazer)
library(glue)
library(frailtypack)

### Loading data
load(file="./sm_lsq_data.rda")

### Main Analyses

## Testing Hypothesis 1. Logistic Regression. 
logistic_model <- glm(RIC_Coalition ~ Ideological_Dispersion + GDP + 
                        Inflation +
                        Pres_Approval + Recess, family="binomial", 
                      data = sm_lsq)
summary(logistic_model)

# Table 2. The Effect of Ideologically-Heterogeneous 
# Cabinets on Policy Monitoring 
# (Dependent Variable: RIC Coalition) 
stargazer(logistic_model,
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001), 
          type="text",
          keep=1)

# Exponential value of the coefficient "Ideological Dispersion"
exp(0.256)

## Testing Hypothesis 2. Survival Analysis.
survival_model <- coxph(Surv(start, stop, survival_month == 2) ~ RIC_Coalition + 
                           RIC_Opposition + GDP + Inflation +
                           Pres_Approval +  Ideological_Dispersion + 
                           Independent_Minister + Government_Size + 
                           First_Cabinet + Cumulative_Replace  + Cumulative_RICs  + 
                           Reeligible + Recess + Electoral_Year, 
                         data = sm_lsq)
summary(survival_model)

# Figure 2. Hazard Ratios: Policy Monitoring Between Coalition Partners 
# (RIC Coalition) and the Survival of Ministers
tiff("./Fig_2.tiff", 
     units="in", width=8, height=6, res=300)
ggforest(survival_model, main="", data=sm_lsq, noDigits = 2)
dev.off()




### Supplementary Material 


## Appendix F. Hypothesis 1
# Table F1. The Effect of Ideological Dispersion on Policy Monitoring
stargazer(logistic_model,
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          type="text")


## Appendix G. Descriptive Statistics
# Table G.1. List of Ministries (Portfolios)
unique(sort(sm_lsq$Portfolio_Name))
sm_lsq %>% 
  arrange(Portfolio_Name) %>% 
  distinct(port_abbrev)

# Table G.2. Variables and Descriptive Statistics
stargazer(sm_lsq[c("survival_month", "RIC_Coalition", "RIC_Opposition", "GDP", "Inflation",
                           "Pres_Approval", "Ideological_Dispersion", "Independent_Minister",
                           "Government_Size",
                           "First_Cabinet",  "Cumulative_Replace", "Cumulative_RICs", "Reeligible",
                           "Recess", "Electoral_Year", "Finance", "Policy", "External",
                           "Protests", "Scandals")], 
          type = "text", digits=2, 
          summary.stat = c("mean", "sd", "min", "max", "n"))

## Appendix H. Proportion of RICs by Administration and Portfolio, 
## and Ministers' Average Survival (in Months)
# Figure H.1. Proportion of Months with RICs Between Coalition 
# Members (RIC Coalition = 1) by Administration
fig_h1 <- sm_lsq %>% 
  ggplot(aes(x=Administration, y=perc_month_ric_coal, fill=as.factor(RIC_Coalition), width=.8)) + 
  geom_bar(stat="identity", colour="black", position="dodge") + 
  scale_fill_grey(start = 0.70, end = 0.30) +
  geom_text(aes(label=perc_month_ric_coal), position=position_dodge(width=0.9), vjust=-0.25) +
  ylab("Percentage of Months") + theme_bw() +
  labs(fill = "RIC by Coalition") +
  theme(legend.position="bottom",
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        axis.title.x = element_text(size = 12))

tiff("./Fig_H1.tiff", 
     units="in", width=8, height=5, res=300)
fig_h1
dev.off()

# Proportion of Months with RICs Between Coalition Members
# (RIC_Coalition = 1) (Left-Plot) and Minister's Average Months
# in Office by Portfolio (Right-Plot) (Per Administration)
plots_ric_coal = list()
plots_mean_port = list()
list_do_call_coal = list()
unique_adm = unique(sm_lsq$Administration)
for(i in unique_adm) {
  plots_ric_coal[[i]] <- sm_lsq %>% 
    filter(Administration == i) %>% 
    group_by(port_abbrev) %>% 
    filter(!duplicated(port_perc_month_ric_coal)) %>% 
    ungroup() %>% 
    dplyr::mutate(port_abbrev = factor(port_abbrev, 
                                       levels=unique(port_abbrev[order(RIC_Coalition, -port_adm_perc_month_ric_coal, port_abbrev)]), 
                                       ordered=TRUE)) %>% 
    ggplot(aes(x=port_abbrev, y=port_adm_perc_month_ric_coal, fill=as.factor(RIC_Coalition))) + 
    geom_bar(stat="identity")+ 
    scale_fill_grey(start = 1, end = 0.30) +
    ylab("Percentage of Months") + xlab("Portfolio") +
    geom_text(aes(label=ifelse(RIC_Coalition==1, port_adm_perc_month_ric_coal, "")), 
              hjust=1, colour="white", size=3.5) + 
    theme_few() + scale_y_continuous(expand = c(0, 0), limits=c(0,101)) +
    coord_flip() + 
    ggtitle(glue("{i} (% Months with RIC Coalition)")) +
    theme(legend.position = "none",
          axis.text.y = element_text(size = 10),
          plot.title = element_text(size = 10)) 
  
  plots_mean_port[[i]] <- sm_lsq %>% 
    filter(Administration == i) %>% 
    group_by(port_abbrev) %>% 
    filter(!duplicated(mean_months_port)) %>% 
    ungroup() %>% 
    dplyr::mutate(port_abbrev = factor(port_abbrev, 
                                       levels=unique(port_abbrev[order(mean_months_port, port_abbrev)]), 
                                       ordered=TRUE)) %>% 
    ggplot(aes(x=port_abbrev, y=mean_months_port)) + 
    geom_bar(stat="identity")+ 
    scale_fill_grey(start = 1, end = 0.30) +
    ylab("Average Survival (in Months)") + xlab("") +
    geom_text(aes(label=mean_months_port), 
              hjust=1, colour="white", size=3.5) + 
    theme_few()  + coord_flip() + 
    ggtitle(glue("{i} ({unique(sm_lsq$adm_month_duration[sm_lsq$Administration==i])} month administration)")) +
    theme(legend.position = "none",
          axis.text.y = element_text(size = 10),
          plot.title = element_text(size = 10)) 
  
  list_do_call_coal[[i]] = grid.arrange(plots_ric_coal[[i]], plots_mean_port[[i]], ncol=2)
}

# Figure H.2. FHC I's Administration (1995-1998)
tiff("./Fig_H2.tiff", 
     units="in", width=8, height=4, res=300)
plot(list_do_call_coal[["FHC I"]])
dev.off()

# Figure H.3. FHC II's Administration (1999-2002)
tiff("./Fig_H3.tiff", 
     units="in", width=8, height=4, res=300)
plot(list_do_call_coal[["FHC II"]])
dev.off()

# Figure H.4. Lula I's Administration (2003-2006)
tiff("./Fig_H4.tiff", 
     units="in", width=8, height=4, res=300)
plot(list_do_call_coal[["Lula I"]])
dev.off()

# Figure H.5. Lula II's Administration (2007-2010)
tiff("./Fig_H5.tiff", 
     units="in", width=8, height=4, res=300)
plot(list_do_call_coal[["Lula II"]])
dev.off()

# Figure H.6. Dilma I's Administration (2011-2014)
tiff("./Fig_H6.tiff", 
     units="in", width=8, height=4, res=300)
plot(list_do_call_coal[["Dilma I"]])
dev.off()

## Appendix I. Survival Analysis
# Table I.1. Hazard Rates: The Effect of Policy Monitoring 
# (RIC Coalition) on the Survival of Individual Ministers
stargazer(survival_model, 
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          type="text")

## Appendix J. Model Fit and Proportional Hazards Assumption Tests
log_negative_log <- coxph(Surv(start, stop, survival_month == 2) ~ 
                                  strata(RIC_Coalition) + GDP + Inflation +
                                  Pres_Approval + RIC_Opposition + 
                                  Ideological_Dispersion + 
                                  Independent_Minister + Government_Size + 
                                  First_Cabinet + Cumulative_Replace  + Cumulative_RICs  + 
                                  Reeligible + Recess + 
                                  Electoral_Year,
                                data = sm_lsq)

# Figure J.1. "Log-negative-log" of the Survival Function
tiff("./Fig_J1.tiff", 
     units="in", width=6, height=5, res=300)
ggsurvplot(survfit(log_negative_log), data = sm_lsq, 
           censor=FALSE,
           fun = "cloglog",
           palette = "npg",
           xlab = "ln(time)",
           ylab = "log[-log(Suurvival Probability)]",
           conf.int = F,
           legend.title ="",
           xlim = c(1,50),
           ggtheme = theme_bw()) 
dev.off()

# Figure J.2.  Schoenfeld Residuals Test 
cox.zph.fit <- cox.zph(survival_model)

tiff("./Fig_J2.tiff", 
     units="in", width=14, height=10, res=300)
ggcoxzph(cox.zph.fit, font.main = 10, font.y=9, font.x=10) 
dev.off()

# Table J.1. Interaction of RIC Coalition with Time Variable "Month"
interaction_model <- coxph(Surv(start, stop, survival_month == 2) ~ 
                                RIC_Coalition*Month + GDP + Inflation +
                                Pres_Approval + RIC_Opposition + 
                                Ideological_Dispersion + 
                                Independent_Minister + Government_Size + 
                                First_Cabinet + Cumulative_Replace  + Cumulative_RICs  + 
                                Reeligible + Recess + 
                                Electoral_Year,
                              data = sm_lsq)

# Table J.1. Column "Coefficient"
stargazer(survival_model,interaction_model, 
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          type="text", 
          keep=c(1,16))

# Table J.1. Column "exp(coef)"
stargazer(survival_model,interaction_model, 
          star.cutoffs = NA,
          type="text", 
          keep=c(1,16),
          apply.coef = exp)

## Appendix K. Robustness Checks
# Table K.1. Benchmark Model (Without Controls)
benchmark_model <- coxph(Surv(start, stop, survival_month == 2) ~ 
                             RIC_Coalition,
                           data = sm_lsq)

stargazer(benchmark_model, 
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          type="text")

# Table K.2. Fixed Effects (FE) Models
month_fe <- coxph(Surv(start, stop, survival_month == 2) ~ RIC_Coalition + 
                                    RIC_Opposition + GDP + Inflation +  Pres_Approval +
                                    Ideological_Dispersion + Independent_Minister + 
                                    Government_Size + First_Cabinet + 
                                    Cumulative_Replace  + 
                                    Cumulative_RICs + Reeligible + 
                                    Recess + Electoral_Year + 
                                    as.factor(Month),
                                  data = sm_lsq)
summary(month_fe)

portfolio_fe <- coxph(Surv(start, stop, survival_month == 2) ~ RIC_Coalition + 
                                   RIC_Opposition + GDP + Inflation +  Pres_Approval +
                                   Ideological_Dispersion + Independent_Minister + 
                                   Government_Size + First_Cabinet + 
                                   Cumulative_Replace  + 
                                   Cumulative_RICs + Reeligible + 
                                   Recess + Electoral_Year +
                                   as.factor(Portfolio_Name),
                                 data = sm_lsq)
summary(portfolio_fe)

coalition_fe <- coxph(Surv(start, stop, survival_month == 2) ~ RIC_Coalition + 
                                   RIC_Opposition + GDP + Inflation +  Pres_Approval +
                                   Ideological_Dispersion + Independent_Minister + 
                                   Government_Size + First_Cabinet + 
                                   Cumulative_Replace  + 
                                   Cumulative_RICs + Reeligible + 
                                   Recess + Electoral_Year +
                                   as.factor(Coalition),
                                 data = sm_lsq)
summary(coalition_fe)

administration_fe <- coxph(Surv(start, stop, survival_month == 2) ~ RIC_Coalition +
                                  RIC_Opposition + GDP + Inflation +  Pres_Approval +
                                  Ideological_Dispersion + Independent_Minister + 
                                  Government_Size + First_Cabinet + 
                                  Cumulative_Replace  + 
                                  Cumulative_RICs + Reeligible + 
                                  Recess + Electoral_Year +
                                  as.factor(Administration),
                                data = sm_lsq)
summary(administration_fe)

stargazer(month_fe, portfolio_fe, coalition_fe, administration_fe,
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          type="text",
          keep=c(1:14))

## Table K.3. Model Extensions: Policy Areas, Protests, and Scandals
rob_areas <- coxph(Surv(start, stop, survival_month == 2) ~ RIC_Coalition + 
                          RIC_Opposition + GDP + Inflation + Pres_Approval +
                          Ideological_Dispersion + Independent_Minister + 
                          Government_Size + First_Cabinet +
                          Cumulative_Replace  + Cumulative_RICs + 
                          Reeligible + Recess + Electoral_Year +
                          Finance + Policy + External, data=sm_lsq)
summary(rob_areas)

rob_prot_scans <- coxph(Surv(start, stop, survival_month == 2) ~ RIC_Coalition + 
                               RIC_Opposition + GDP + Inflation + Pres_Approval +
                               Ideological_Dispersion + Independent_Minister + 
                               Government_Size + First_Cabinet +
                               Cumulative_Replace  + Cumulative_RICs + 
                               Reeligible + Recess + Electoral_Year +
                               Protests + Scandals, data=sm_lsq)
summary(rob_prot_scans)


stargazer(rob_areas, rob_prot_scans, 
          star.char = c("*", "**", "***"),
          star.cutoffs = c(.05, .01, .001),
          type="text")

## Appendix L. Shared Gamma Frailty Model
# Comparing Results from the Original Cox Model and a 
# Shared Gamma Frailty Model Using Penalized Likelihood 
# Estimation (n.knots = 20, kappa = 1000000)
frailty_pen_model <- frailtyPenal(Surv(start, stop, survival_month == 2) ~ 
                              cluster(Portfolio_Name) + RIC_Coalition + RIC_Opposition + 
                              GDP + Inflation + Pres_Approval + Ideological_Dispersion + 
                              Independent_Minister + Government_Size + 
                              First_Cabinet + Cumulative_Replace  + Cumulative_RICs  + 
                              Reeligible + Recess + 
                              Electoral_Year,n.knots=20, kappa=10000,
                            data = sm_lsq)
frailty_pen_model


